In [1]:
%matplotlib inline
import pandas as pd
import zipfile
import shutil
import urllib2
from urllib2 import urlopen
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.dates as dates
import matplotlib.ticker as tick
import scipy.stats as sp
import statsmodels.api as sm
from pandas.stats.api import ols
from datetime import datetime
from bs4 import BeautifulSoup
from pylab import rcParams
import platform
rcParams['figure.figsize'] = 15, 10
import re
import os
import sys
import glob
import urllib
import HTMLParser
from cStringIO import StringIO

import pyproj
from pyproj import Proj

import gzip

import ftplib
import calendar
import datetime
from datetime import date

import pymodis

In [2]:
import arcpy
arcpy.CheckOutExtension("spatial")
from arcpy import env 
from arcpy.sa import *

In [3]:
print("Operating System " + platform.system() + " " + platform.release())
print("Python Version " + str(sys.version))
print("Pandas Version " + str(pd.__version__))
print("Numpy Version " + str(np.__version__))


Operating System Windows 7
Python Version 2.7.10 (default, May 23 2015, 09:44:00) [MSC v.1500 64 bit (AMD64)]
Pandas Version 0.19.2
Numpy Version 1.12.0

In [4]:
import wellapplication as wa
import UBM

In [5]:
Zonal_HUCS = "H:/GIS/BCM/ZonalExtra.gdb/Zonal_HUC"
Zone_field = "HUC_12"
z_Name = "H:/GIS/BCM/ZonalExtra.gdb"
arcpy.env.overwriteOutput = True

In [ ]:
indata = "H:/GIS/BCM/AvailableWater.gdb"
UBM.zone_gdb(indata, z_Name, Zonal_HUCS, Zone_field)


z_AVWT200401
z_AVWT200402
z_AVWT200403
z_AVWT200404
z_AVWT200405
z_AVWT200406
z_AVWT200407
z_AVWT200408
z_AVWT200409
z_AVWT200410
z_AVWT200411
z_AVWT200412
z_AVWT200501
z_AVWT200502
z_AVWT200503
z_AVWT200504
z_AVWT200505
z_AVWT200506

In [ ]:


In [ ]:


In [ ]:


In [ ]:


In [ ]:
def replace_hdr_file(hdrfile):
    """
    Replace the .hdr file for a .bil raster with the correct data for Arc processing
    Required: hdrfile -- filepath for .hdr file to replace/create
    Output:   None
    """
    # hdr file replacment string
    HDRFILE_STRING = "byteorder M\n\
    layout bil\n\
    nbands 1\n\
    nbits 16\n\
    ncols 6935\n\
    nrows 3351\n\
    ulxmap -124.729583333331703\n\
    ulymap 52.871249516804028\n\
    xdim 0.00833333333\n\
    ydim 0.00833333333\n"
    with open(hdrfile, 'w') as o:
        o.write(HDRFILE_STRING)

In [ ]:
path= "H:/GIS/MODIS2/MODIS.gdb/"
arcpy.env.workspace = path

for rast in arcpy.ListRasters()[:3]:
        dsc = arcpy.Describe(rast)
        print(dsc.baseName)

In [ ]:
for hdrfile in glob.glob("H:/GIS/SNODAS/SNODASUNZIPPED/SNODASUNGZ/*.Hdr"):
    replace_hdr_file(hdrfile)

In [ ]:
sr = arcpy.SpatialReference(4326) #Spatial Projection WGS84

# Set output coordinate system
outCS = arcpy.SpatialReference('NAD 1983 UTM Zone 12N')  #NAD83 Zone 12 Code is 26912

# Set local variables
inMaskData = "H:/GIS/SNODAS/SNODAS.gdb/UT_HUC_area"

# from Table 4 of http://nsidc.org/data/docs/noaa/g02158_snodas_snow_cover_model/
prodcodelong = {1025: 'Precipitation', 1034: 'Snow water equivalent', 1036: 'Snow depth', 
            1038: 'Snow pack average temperature', 1039: 'Blowing snow sublimation',
            1044: 'Snow melt', 1050: 'Snow pack sublimation'}
prodcode = {1025: 'PREC', 1034: 'SWEQ', 1036: 'SNOD', 1038: 'SPAT', 1039: 'BSSB', 1044: 'SNML', 1050: 'SPSB'}

# iterate through data files creating rasters and defining projections
for dtfile in glob.glob("H:/GIS/SNODAS/SNODASUNZIPPED/SNODASUNGZ/*.dat"):
    indata = arcpy.Raster(dtfile)
    arcpy.DefineProjection_management(indata, sr)
    typAbb = prodcode[int(os.path.basename(dtfile)[8:12])]
    
    #indata.save("H:/GIS/SNODAS/SNODAS.gdb/" + typAbb + dtfile[-20:-11]) #YYYYMMDD
    
    # Execute ExtractByMask to clip snodas data to Utah watersheds
    outExtractByMask = arcpy.sa.ExtractByMask(indata, inMaskData)

    # Determine the new output feature class path and name = productcode + YYYYMMDD and save to a geodatabase
    outrast = "H:/GIS/SNODAS/SNODAS.gdb/" + typAbb + dtfile[-20:-11] #os.path.join(outWorkspace, rast)
    
    # Print name to veryfy save
    print(typAbb + dtfile[-20:-11])
    
    # Project Raster to UTM Zone 12
    arcpy.ProjectRaster_management(outExtractByMask, outrast, outCS, 'BILINEAR', '1000',\
                                   'WGS_1984_(ITRF00)_To_NAD_1983', '#', '#')

From http://nsidc.org/data/docs/noaa/g02158_snodas_snow_cover_model/, the file abbreviations are as follows:

  • `RAIN` = `Wet Precip`
  • `SWEQ` = `Snow Water Equivalent`
  • `SNOD` = `Snow Depth`
  • `SPAT` = `Snow Pack Average Temp`
  • `BSSB` = `Blowing Snow Sublimation`
  • `SNML` = `Snowmelt`
  • `SPSB` = `Snow Pack Sublimation`

In [ ]:
prodcode = {'us_ssmv11038wS__A':'SPAT', 'us_ssmv11044bS__T':'SNML', 'us_ssmv11050lL00T':'SPSB', 
            'us_ssmv11034tS__T':'SWEQ', 'us_ssmv01025SlL00':'RAIN', 'us_ssmv01025SlL01':'SNOW',
            'us_ssmv11036tS__T':'SNOD', 'us_ssmv11039lL00T':'BSSB'}

path = "H:/GIS/SNODAS/geotifSNODAS/SNDS/"
for filename in os.listdir(path):
    if filename.startswith("us_ssmv"):
        code = prodcode[filename[0:17]]
        yrsrt = filename.find('TNATS') + 5
        yr = filename[yrsrt:yrsrt+4]
        mo = filename[yrsrt+4:yrsrt+6]
        dy = filename[yrsrt+6:yrsrt+8]
        os.rename(os.path.join(path, filename), os.path.join(path,code+yr+mo+dy+filename[-4:]))

In [ ]:
def mergeRasts(path, data_type = 'AET', monthRange = [1,4], yearRange = [2000,2001]):
    arcpy.env.workspace = path
    print(arcpy.ListRasters())
    
    for y in range(yearRange[0],yearRange[-1]+1): #set years converted here
        for m in range(monthRange[0],monthRange[-1]+1): #set months converted here   
            nm = data_type + str(y) + str(m).zfill(2)
            rlist=[]
            for rast in arcpy.ListRasters(nm+'*'): 
                rlist.append(rast)
            print(rlist)

In [ ]:
path="H:/GIS/MODIS2/MODIS.gdb"
mergeRasts(path)

In [ ]:
path = "H:/GIS/SNODAS/geotifSNODAS/SNDS/"
g = {}
from arcpy.sa import *
arcpy.env.workspace = path
arcpy.env.overwriteOutput = True

code = 'PREC'

for y in range(2003,2016):
    for m in range(1,13):
        g[code+str(y)+str(m).zfill(2)] = []
        for name in sorted(glob.glob(path+code+'*.tif')):
            rast = os.path.basename(name)
            if rast[0:4] == code and int(rast[4:8]) == y and int(rast[8:10]) == m:
                g[code+str(y)+str(m).zfill(2)].append(rast)
            else:
                pass
        if len(g[code+str(y)+str(m).zfill(2)])>0:
            print(g[code+str(y)+str(m).zfill(2)])
            calc = CellStatistics(g[code+str(y)+str(m).zfill(2)], statistics_type = "SUM", ignore_nodata="DATA")
            calc.save("H:/GIS/SNODAS/SNODAS.gdb/"+rast[0:4]+str(y).zfill(2)+str(m).zfill(2)+"SUM")

In [ ]:


In [ ]:
sr = arcpy.SpatialReference(4326) #Spatial Projection WGS84

# Set output coordinate system
outCS = arcpy.SpatialReference('NAD 1983 UTM Zone 12N')  #NAD83 Zone 12 Code is 26912

# Set local variables
inMaskData = "H:/GIS/SNODAS/SNODAS.gdb/UT_HUC_area"

# from Table 4 of http://nsidc.org/data/docs/noaa/g02158_snodas_snow_cover_model/
prodcodelong = {1025: 'Precipitation', 1034: 'Snow water equivalent', 1036: 'Snow depth', 
            1038: 'Snow pack average temperature', 1039: 'Blowing snow sublimation',
            1044: 'Snow melt', 1050: 'Snow pack sublimation'}
prodcode = {1025: 'PREC', 1034: 'SWEQ', 1036: 'SNOD', 1038: 'SPAT', 1039: 'BSSB', 1044: 'SNML', 1050: 'SPSB'}

# iterate through data files creating rasters and defining projections
for dtfile in glob.glob("H:/GIS/SNODAS/SNODASUNZIPPED/SNODASUNGZ/*.dat"):
    indata = arcpy.Raster(dtfile)
    arcpy.DefineProjection_management(indata, sr)
    typAbb = prodcode[int(os.path.basename(dtfile)[8:12])]
    
    # Determine the new output feature class path and name = productcode + YYYYMMDD and save to a geodatabase
    outrast = "H:/GIS/SNODAS/SNODAS.gdb/" + typAbb + dtfile[-20:-11] #os.path.join(outWorkspace, rast)
    
    # Print name to veryfy save
    print(typAbb + dtfile[-20:-11])
    
    # Project Raster to UTM Zone 12
    arcpy.ProjectRaster_management(outExtractByMask, outrast, outCS, 'BILINEAR', '1000',\
                                   'WGS_1984_(ITRF00)_To_NAD_1983', '#', '#')

In [ ]:
arcpy.env.overwriteOutput = True
path="H:/GIS/SNODAS/SNODASproj.gdb/"
arcpy.env.workspace = path
for name in arcpy.ListRasters():

In [ ]:
monthRange = [1,12] 
yearRange = [2003,2016]

g = {}
path="H:/GIS/SNODAS/SNODASproj.gdb/"
arcpy.env.workspace = path
arcpy.env.overwriteOutput = True


for y in range(yearRange[0],yearRange[1]+1): #set years converted here
    for m in range(monthRange[0],monthRange[1]+1): #set months converted here
        my = str(y)+str(m).zfill(2)
        newdn = 'TPPT' + my
        try:
            calc = Plus('SNOW'+ my +'SUM', 'RAIN'+ my +'SUM')
            calc.save(newdn+'SUM')
            print(newdn)
        except(RuntimeError):
            pass

In [ ]:
from arcpy.sa import *

monthRange = [1,12] 
yearRange = [2005,2005]

path = "H:/GIS/Calc.gdb/"
path1 = "H:/GIS/SNODAS/SNODASproj.gdb/"
path2 = "H:/GIS/MODIS/MODOUT.gdb/"
arcpy.env.workspace = path
arcpy.env.overwriteOutput = True
area = 'H:/GIS/NHD_UT_Proj.gdb/UT_HUC_area'
arcpy.env.mask = area

for y in range(yearRange[0],yearRange[1]+1): #set years converted here
    for m in range(monthRange[0],monthRange[1]+1): #set months converted here
        my = str(y)+str(m).zfill(2)
        newdn = 'AVLW' + my
        rain = (path1 + 'RAIN'+ my +'SUM')
        snowMelt = (path1 + 'SNML'+ my +'SUM')
        actEvap = (path2 + 'AET'+ my)
        
        avail = (Con(IsNull(rain),0, rain)) + (Con(IsNull(snowMelt),0, snowMelt)) - (Con(IsNull(actEvap),0, actEvap))
        avail = (available < 0, 0, available)
        available.save(newdn)
        print(newdn)

In [ ]: